function  [uc vc] = Find_StablePoint ( a, b, umax, umin )
% find the equilibrium point

tmax=1000000; % the number of calculation
uc=(umax+umin)/2;
fc=a/(1+uc^b)-(a/uc-1)^(1/b);

for i=1:1:tmax;
    if fc>0
       umax=(umax+umin)/2;
       uc=(umax+umin)/2;
       fc=a/(1+uc^b)-(a/uc-1)^(1/b);
    elseif fc<0
       umin=(umax+umin)/2;
       uc=(umax+umin)/2;
       fc=a/(1+uc^b)-(a/uc-1)^(1/b);
    else
    break;
    end
end
vc=a/(1+uc^b);
hold on;
plot(uc,vc,'ro');

end

